/*************************************************************/
/*    Source for the stress of a neo-hookean finite strain   */
/*    Initially based on ComputeFiniteStrainElasticStress    */
/*                  Status: not working                      */
/*                                                           */
/* Copyright (C) 2018  eDgar                                 */
/*                                                           */
/* This program is free software: you can redistribute it    */
/* and/or modify it under the terms of the GNU General       */
/* Public License as published by the Free Software          */
/* Foundation, either version 3 of the License, or (at your  */
/* option) any later version.                                */
/*                                                           */
/* This program is distributed in the hope that it will be   */
/* useful, but WITHOUT ANY WARRANTY; without even the        */
/* implied warranty of MERCHANTABILITY or FITNESS FOR A      */
/* PARTICULAR PURPOSE.  See the GNU General Public License   */
/* for more details.                                         */
/*                                                           */
/* You should have received a copy of the GNU General Public */
/* License along with this program.  If not, see             */
/* <http://www.gnu.org/licenses/>.                           */
/*************************************************************/

#include "ComputeStressNeohookean.h"

registerMooseObject("spiderApp", ComputeStressNeohookean);


// Just an idea to stop copying the code from ComputeElasticityTensorNeohookean
// template <>
// InputParameters
// validParams<ElementAverageValue>()
// {
//   InputParameters params = validParams<ElementIntegralVariablePostprocessor>();
//   return params;
// }


template <>
InputParameters
validParams<ComputeStressNeohookean>()
{
  InputParameters params = validParams<ComputeStressBase>();
  params.addClassDescription
    ("Compute neo-Hookean stress using elasticity for "
     "finite strains");
  /*
   * This code is commented, because I am going to hard code
   * parameters to see if that fixes issues with the coupon
   *
   * ///////////////////////////////////////
   * // Repeated code from                //
   * // ComputeElasticityTensorNeohookean //
   * ///////////////////////////////////////
   *
   * // Give a default value to \mu (set -1; if the value is
   * // not set to a positive value by the user, there is a
   * // check below which throws an error)
   * params.addParam<Real>
   *   ("mu",
   *    -1.0,
   *    "Lamé second constant (shear) of the material.");
   * // Allow lambda to be an input parameter
   * params.addParam<Real>
   *   ("lambda",
   *    "Lame's first constant (non-shear) of the material.");
   *
   * // Allow Poisson's ratio and Young's modulus as input
   * // parameters
   * params.addParam<Real>
   *   ("nu", "Poisson's ratio of the material.");
   * params.addParam<Real>
   *   ("E", "Young's modulus of the material.");
   */

  return params;

  ///////////////////////////////////////
  // -- end Repeated code              //
  ///////////////////////////////////////
}

ComputeStressNeohookean::ComputeStressNeohookean(
    const InputParameters & parameters)
  : ComputeStressBase(parameters),
    GuaranteeConsumer(this),

    /*
     * This code is commented, because I am going to hard code
     * parameters to see if that fixes issues with the coupon
     *
     * ///////////////////////////////////////
     * // Repeated code from                //
     * // ComputeElasticityTensorNeohookean //
     * ///////////////////////////////////////
     *
     * // Get whether or not the parameters were given
     * _lambda_set(parameters.isParamValid("lambda")),
     * _mu_set(parameters.isParamValid("mu")),
     * _nu_set(parameters.isParamValid("nu")),
     * _E_set(parameters.isParamValid("E")),
     *
     * // Take the Lamé constants from
     * // the GetPot (.i) file
     * _lambda(_lambda_set ? getParam<Real>("lambda") : -1),
     * // The same for `mu'
     * _mu(_mu_set ? getParam<Real>("mu") : -1),
     *
     * // `nu' and `E' are parameters, because I don't need to
     * // store them for anything else (so far).
     * _nu(_nu_set ? getParam<Real>("nu") : -1),
     * _E(_E_set ? getParam<Real>("E") : -1),
     */
     // Take the deformation gradient from
     // ComputeFiniteStrain
     _deformation_gradient(getMaterialProperty<RankTwoTensor>
     ("deformation_gradient"))
    /*
     * ///////////////////////////////////////
     * // -- end Repeated code              //
     * ///////////////////////////////////////
     *
     * // "You can declare material property for J_1 in the
     * // ComputeElasticityTensorNeohookean.C and get material
     * // property J_1 in ComputeStressNeohookean.C." --Wen Jiang
     * //
     * // I hope that this is the way to use a property that I
     * // declared in ComputeElasticityTensorNeohookean
     * _J_1(getMaterialProperty<Real>
     *      ("jacobian_inverse")),
     * _logJ(getMaterialProperty<Real>
     *       ("jacobian_log")),
     *
     * // See note about "Rotate the stress" in source code
     * // (following ComputeFiniteStrainElasticStress.C and
     * // ComputeFiniteStrain::computeQpIncrements).
     * _rotation_increment(getMaterialProperty<RankTwoTensor>
     *                     ("rotation_increment"))
     */
{
  /*
   * ///////////////////////////////////////
   * // Repeated code from                //
   * // ComputeElasticityTensorNeohookean //
   * ///////////////////////////////////////
   * // Neo-hookean material model
   *
   * // Now we are going to initialise variables for the analysis
   * // TODO:
   *
   * // Make sure that only one of these are given:
   * // 1. Young's modulus and the Poisson's ratio
   * // 2. Lamé's constants lambda and mu
   * if (! ( (_E_set && _nu_set) || (_lambda_set && _mu_set) ) )
   *   mooseError
   *     ("One of these combinations is needed:\n"
   *      "(1) Young's modulus (E) and the Poisson's ratio (nu)\n"
   *      "(2) Lamé's constants lambda and mu\n"
   *      "in the input file for the neo-Hoookean\n"
   *      "'" + name() + "'.");
   * else if (_E_set && _nu_set) {
   *     // If both Young's modulus and Poisson's ratio are
   *     // given, use them to calculate Lamé's constants
   *     if (_nu_set && (_nu <= -1.0 || _nu >= 0.5))
   *       mooseError
   *         ("Poissons ratio must be greater than -1 and less"
   *          " than 0.5 in material '" + name() + "'.");
   *     if (_E_set && _E <= 0.0)
   *       mooseError
   *         ("Youngs modulus must be positive in material '"
   *          + name() + "'.");
   *
   *     // Value of Lamé constants in terms of Young's modulus
   *     // and Poisson's ratio
   *     _mu = _E / (2 * (_nu + 1));
   *     // // 3 multiplications, 1 division, 1 addition, 1 subtraction:
   *     // _lambda = E * _nu/ ((_nu + 1) * (1 - 2 * _nu));
   *     // 3 multiplications, 1 division, 1 subtraction, less clear:
   *     _lambda = 2 * _mu * _nu / (1 - 2 * _nu);
   *   }
   * else if (_lambda_set && _mu_set) {
   *   // Make sure that the values of Lamé's constants are
   *   // not both negative or zero.
   *   // Holzapfel, G. A., Nonlinear solid mechanics : a
   *   // continuum approach for engineering (2000), :
   *   // Wiley. p.250
   *   if ((_lambda < 0.0 || _mu < 0.0) ||
   *       (_lambda == 0.0 && _mu == 0.0))
   *     mooseError("The Lamé constants must be positive in material '"
   *                + name() + "'.");
   * }
   * ///////////////////////////////////////
   * // -- end Repeated code              //
   * ///////////////////////////////////////
   */
}

void
ComputeStressNeohookean::initialSetup()
{
  if (!hasGuaranteedMaterialProperty
      (_elasticity_tensor_name, Guarantee::ISOTROPIC))
    mooseError("ComputeStressNeohookean can only be used "
               "with elasticity tensor materials "
               "that guarantee isotropic tensors.");
}

void
ComputeStressNeohookean::initQpStatefulProperties()
{
  ComputeStressBase::initQpStatefulProperties();
}

void
ComputeStressNeohookean::computeQpStress()
{
  // \LaTeX:
  /** \begin{align}
   *    \sigma_{ij} = \frac{1}{J}\left(
   *      \lambda_{0}\cdot\ln{\left(J\right)}\cdot\delta_{ij}
   *       + \mu_{0}\cdot \left(B_{ij}
   *          - \delta_{ij}\right)
   *    \right)
   *  \end{align}
   */

  ////////////////////////////
  // Mathematical constants //
  ////////////////////////////

  // Identity matrix
  // 1  0  0
  // 0  1  0
  // 0  0  1
  // initIdentity() in
  // [[moose/framework/src/utils/RankTwoTensor.C]]
  const RankTwoTensor dij(RankTwoTensor::initIdentity);

  /////////////////////////
  // Material parameters //
  /////////////////////////

  // Hard-code Young's modulus for testing
  const Real E = 30e6;
  // Hard-code Poisson's ratio for testing
  const Real nu = 0.4;

  // Value of Lamé constants in terms of Young's modulus
  // and Poisson's ratio
  const Real mu = E / (2.0 * (nu + 1.0));
  // // 3 multiplications, 1 division, 1 subtraction, less clear:
  // _lambda = 2.0 * _mu * _nu / (1.0 - 2.0 * _nu);
  // 3 multiplications, 1 division, 1 addition, 1 subtraction:
  const Real lambda = E * nu/ ((nu + 1) * (1 - 2 * nu));

  /////////////////////////
  // Kinematic variables //
  /////////////////////////

  // Deformation gradient tensor F = 1 + du_{i}/dX_{j}
  // Just copy the value for the given quadrature point
  const RankTwoTensor F = _deformation_gradient[_qp];
  // Jacobian (determinant of the deformation gradient)
  // Instead of this:
  // Real _J = _deformation_gradient[_qp].det();
  // Do this:
  const Real J = F.det();
  // Inverse of the deformation Jacobian
  const Real J_1 = 1 / J;
  // Natural logarithm of the deformation Jacobian
  const Real logJ = log(J);

  // Left Green-Cauchy deformation tensor
  const RankTwoTensor Bij = F * F.transpose();

  _stress[_qp] = J_1 * (lambda * logJ * dij + mu * (Bij - dij));

  /*
  * // Rotate the stress state to the current configuration
  * // (following ComputeFiniteStrainElasticStress.C and
  * // ComputeFiniteStrain::computeQpIncrements).
  * // We assume that the calculations of stress are done
  * // without the rotation portion of F = R·U, and _stress
  * // actually needs it. It would work awfully without it for
  * // linear behaviour and not work at all for non-linear
  * // materials.
  * _stress[_qp] = _rotation_increment[_qp] \
  *   * _stress[_qp] * _rotation_increment[_qp].transpose();
  * // _stress[_qp] = RankTwoTensor(RankTwoTensor::initIdentity) * 2;
  */

  // "Jacobian_mult is a fourth order tensor that is the
  // partial derivative of your stress with respect to the
  // strain" --Michael Tonks
  //
  // The question was answered by Dr. Lujan with the help of
  // [[cite:BonetWood99][Eq. 5.13]]: \[\dot{\sigma_{ij}} =
  // c_{ijkl} : \dot{d_{kl}}\], where \(d_{kl}\) is the
  // push-forward of \(\dot{E}\), which implies that
  // \(c_{ijkl} = \frac{\dot{\sigma_{ij}}}{\dot{\d_{kl}}}\)
  //
  // Compute dstress_dstrain
  _Jacobian_mult[_qp] = _elasticity_tensor[_qp];
  // Also see: FiniteStrainHyperElasticViscoPlastic.C
  // _Jacobian_mult[_qp] = _tan_mod * _dfe_df * _df_dstretch_inc;

  // // WAIT: testing
  // // From ComputeFiniteStrainElasticStress.C, (and defined
  // // in ComputeFiniteStrain.C)
  // // Assign value for elastic strain, which is equal to the
  // // mechanical strain
  // _elastic_strain[_qp] = _mechanical_strain[_qp];
}
